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ABSTRACT 

We examine the velocity distribution function (VDF) in dark matter halos from Milky Way to 
cluster mass scales. We identify an empirical model for the VDF with a wider peak and a steeper 
tail than a Maxwell-Boltzmann distribution, and discuss physical explanations. We quantify sources 
of scatter in the VDF of cosmological halos and their implication for direct detection of dark matter. 
Given modern simulations and observations, we find that the most significant uncertainty in the VDF 
of the Milky Way arises from the unknown radial position of the solar system relative to the dark 
matter halo scale radius. 

Subject headings: dark matter — galaxy: halo 
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1. INTRODUCTION 

Dark matter is the dominant component of mat- 
ter in the Universe, and the key to the formation of 
large-scale and galactic structures. Modern cosmo- 
logical observations suggest that dark matter is com- 
posed of a yet- unidentified elementary particle (e.g. lFengl 
2010). However, direct evidence for dark matter par- 
ticles has proved elusive. Experiments that search 
for Weakly Interacting Massive Particles (WIMPs), one 
of the most plausible particle dark matter candidates, 
seek to identify the scattering of a WIMP with a nu 
cleus in an underground detector (Bernabei et al.1 12008 
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and eventually measuring, the WIMP mass and cross sec- 
tion requires a precise understanding of the dark matter 
spatial and veloci ty distribution at the Earth's location 
in the Milky Way dStrigari fc Trottal l200ll IMcCabe1[20Tol: 
IReed et al.ll2T)Tll : IGreenll2012h . 

Dark matter is distributed in halos extending beyond 
the visible components of galaxies; many statistical prop- 
erties including the formation and structure of these ha- 
los have been well characterized by simulations. Despite 
the diversity in the merger and accretion histories of dark 
matter halos of different masses, cosmological simula- 
tions have long sugg ested near universality in the den- 
sity profiles of halos (jNavarro et al.lll996l 119971) . There 
have been several attempts to connect this universality in 
the density profile to t he dark matter Velocity Distribu- 
tion Function (VDF) (lHansen et al J 12003 iKuhlen et all 
12010: INavarro et all 12010) . However, there is no well- 
established model or description for the VDF that has 
been rigorously tested with cosmological simulations. 

Both the implications for direct detection and the quest 
for a theoretical understanding of the phase-space dis- 
tribution in dark matter halos motivate a study of the 
VDF. Under specific, and perhaps too stringent, assump- 
tions, including isolation, equilibrium, spherical symme- 
try, and isotropy, the VDF may be determined uniquely 
from the density profile. For example, with all assump- 
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tions named above and a known density profile, the er- 
godic distribution function can be calculated using Ed- 
dington's formula (Eddington 1916). Although useful as 
an analytic framework, these assumptions are unlikely to 
strictly hold for halos formed via hierarchical merging. 

In absence of an understanding from first principles, 
a practical approach to study the VDF involves appeal- 
ing directly to dark matter halos with a wide range of 
physical properties in cosmological simulations. Quan- 
tifying the VDF directly from cosmological simulations 
would provide a better empirically-motivated framework 
to predict signals in direct detection experiments. Fur- 
thermore, with a parametrized VDF, it becomes more 
tractable to study the relations between the VDF and 
other physical quantities of the halos, such as mass, den- 
sity profile, shape, and formation history. 

In this study, we use a suite of dark matter halos from 
cosmological simulations to study the VDFs at differ- 
ent radii of these halos. We identify a similarity in 
VDFs among a wide range of halos with different masses, 
concentrations, and other physical quantities, that de- 
pends primarily on r/r s , the radius at which it is mea- 
sured divided by the scale radius of the density pro- 
file. We further noti ce that neither standa rd Maxwell- 
Boltzmann models (|Lewin fc Smith! [1996) nor models 
that have been previously prop osed to describe collision- 
less structures (lHansen et all 120061: IKuhlen et all 120101 : 
INavarro et al.ll2010l ) are able to provide an adequate de- 
scription of cosmological VDFs. Instead, we describe the 
distribution of the norm of velocity (in the Galactic rest 
frame) more accurately with an empirical model: 



/(|v|) 



Aexp(-\v\/v ) (y, 
0, 



2 

esc 



, < |v| < w osc 
otherwise, 

where the normalization constant A is chosen such that 
the integral 4ir Jq BBC v 2 f(v)dv equals the number of par- 
ticles in the region of interest. Note that in this pa- 
rameterization the VDF approaches an exponential dis- 
tribution instead of a Gaussian distribution at the low- 
velocity end. With this model, we quantify the scatter 
in the VDF from a variety of sources, including halo-to- 
halo scatter, scatter from finite particle sampling, and 
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scatter from the uncertain position of the Earth within a 
given halo. We further identify the largest uncertainties 
that currently exist in our understanding of the VDF 
at the location of the Earth in our Galaxy, and quan- 
tify their relevance for inferences from direct detection 
experiments. 

2. UNIVERSAL VELOCITY DISTRIBUTION IN 
SIMULATIONS 

To identify the relevant physical quantities which af- 
fect the VDF and to quantify scatter in the distribu- 
tions among different halos in cosmological simulations, 
we must examine a large number of halos across a wide 
range of mass. We also need high resolution to reduce 
sampling error and distinguish differences in VDFs for 
different parameters. 

In this study, we use halos from the Rhapsody and 
Bolshoi simulations; state-of-the-art dark-matter-only 
simulations with high mass resolution. Rhapsody con- 
sists of re-simulations of 96 massive cluster-size halos 
with M vir = 10 14 - 8±0 05 M Q /i" 1 . The particle mass is 
1.3 x 10 8 M Q h~ 1 , resulting in ~ 5 x 10 6 particles in each 
halo. This simulation set currently comprises the largest 
number of halos simulated with this many particles in a 
narrow mass bin (Fig. 1 of lWu et aLll2012[ ). Bolshoi is 
a full cosmological simulation, with similar mass resolu- 
tion, 1.3 x lO 8 M0ft _1 . For detailed descripti ons of the 
Rhap sody and Bolsho i simu lations, refer to IWu et "all 
(2012) and lKlvpin et all (|2011f ) respectively. 

We use the phase-spa ce halo finder Rock- 
STAR (jBehroozi et al.1 1201 lh to identify host halos 
at z = 0. The masses and radii of the halos are 
defined by the spherical overdensity of virialization, 
M{< r vir ) = ^rjJjj.AvirPc, where A vlr = 94 and p c is 
the critical density. We examine the VDFs at a range 
of radii. A VDF at radius r uses all particles within a 
spherical shell centered at the halo center with the inner 
and outer radii of 10 °' 05 r, so that the ratio of the shell 
width to the radius is fixed. In each shell, we assign 
the escape velocity (v cac ) as the spherically-averaged 
v csc of all particles in the shell. We have verified 
that f csc determined from this method is consistent 
with the same quantity deduced from the best-fitting 
spherically-averaged smooth density profile. 

We fit each halo with an NFW density profile, 

(r/r s )(l + r/r s y 

where r s is the scale radius at which the log-log slope is 
—2. The fit uses maximum- likelihood estimation based 
on particles within r V i r . The halo concentration is defined 
as c = r v j r / v s . 

Fig. [1] shows the VDF at different values of r/r s . The 
value of r/r s affects the shape of VDF dramatically. The 
peak of the distribution is a strong function of r/r s . If 
instead the velocity is normalized by the circular veloc- 
ity at each radius rather than the escape velocity, this 
trend will be slightly weakened but still significant. This 
trend in r/r s is not surprising because the VDF heav- 
ily depends on the gravitational potential. If the density 
profiles of simulated halos can be described by the NFW 
profile, which is a function of r/r s only (up to a nor- 
malization constant), the VDF should mostly depend on 
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Figure 1. Solid colored lines show the stacked velocity distribu- 
tion for 96 halos in RHAPSODY, at different values of r/r s : (from left 
to right) 0.15 (blue), 0.3 (red), 0.6 (green), 1.2 (magenta). Bands 
show the 68% halo-to-halo scatter in those VDFs. Dashed and dot- 
ted colored lines indicate the same values of r/r a in BOLSHOI with 
halos of Mvir ~ 10 12 and W ri Mq^ 1 respectively. The VDFs of 
low-mass halos are cut at the head and tail due to limited particle 
number, and their scatter is not shown. The SHM (no = 220 km/s 
and v CBC = 544 km/s) is shown for comparison (black). 

r/r s until the isolated NFW potential breaks down at 
large radius. 

The above trend is robust for halo masses down to 
~ 10 12 M Q , as shown by the Bolshoi simulation in 
Fig. Q] The scatter of the VDFs in the low-mass halos 
considered is somewhat larger due to resolution. How- 
ever, when the high-mass halos are downsampled to have 
the same particle number, the spreads in the stacked 
VDF are comparable to the low-mass halos. We further 
investigated the impact of a variety of parameters char- 
acterizing the halo on the shape of the VDF, and found 
that for a fixed value of r/r s , the halo- to- halo scatter in 
the VDFs is not significantly reduced when binning on 
concentration, shape, or formation history. A detailed 
discussion on this halo-to- halo scatter is in Section 2) 

3. MODELS OF THE VELOCITY DISTRIBUTION 
FUNCTION 

The dark matter velocity distribution in halos is set 
by a sequence of mergers and accretion. The pro- 
cess of violent relaxation (|Lvnden-Belll I1967T) may be 
responsible for the resulting near-equilibrium distribu- 
tions observed in dark matter halos and in galaxies. 
These near-equilibrium distr ibutions explain why ex- 
isting VDF models (see e.g. iFrandsen et al.l [20121 ). in- 
cluding the Standard Halo Model (SHM), King model, 
the double power-law model, and the Tsallis model, 
are all variants of the Maxwell-Boltzmann distribu- 
tion. Recent studies have shown that the widely- 
used SHM, which is a Maxwell-Boltzmann distribu- 
tion with a cut-off put in by hand, is inconsistent 
with th e VDF found in a handful of individual simu- 
lations dStiff fc Widrowl [200l iVpgelsberger et ail [200l 
IKuhlen et alJl2010HPurcell et al.ll2012l) and in the study 
of rotation curve data (jBhattachariee et ai] |2012). The 
double power-law model was proposed to suppress the 
tail of the distribu tion, by raising th e SHM to the power 
of a parameter k (|Lisanti et al.ll20ill) . The Tsallis model 
replaces the Gaussian in Maxwell-Boltzmann distribu- 
tion with a g-Gaussian, which approaches to a Gaussian 
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Figure 2. The VDF for one representative dark matter halo 
in Rhapsody (histogram), along with the best fits using Eq. JT} 
with (vo/vescP) = (0.13,0.78) (black, x 2 = 0.59), SHM (blue, 
9.67), the double power-law model (cyan, 9.47), the Tsallis model 
(green, 1.99), and the analytic VDFs from Eddington's formula 
with isotropic assumption (red dash, 8.48), Osipkov-Merritt (ma- 
genta dash, 6.41), and constant = 1/2 (yellow dash, 11.8). The 
y-axis is in log scale in the main figure and linear in the inset. 

as q — > 1 ( Vcrga dos et al.l |2008|) . It was argued that 
the Tsall is model provide s better fit to simulations with 
baryons (|Ling et al.ll2010t ) , although this conclusion may 
be affected by the relatively low resolution of the simu- 
lations. 

In contrast, our empirical model, Eq. ([1]), is not based 
on a Gaussian distribution but rather on an exponential 
distribution. It also has a power-law cut-off in (binding) 
energy. Fig. [2] shows the VDF in a simulated halo, along 
with the best fit from Eq. ([T]) and the best fits from other 
conventional models. All the best-fit parameters are ob- 
tained from the maximum-likelihood estimation in the 
range of (O^esc)- The fits using Eq. ([1]) are statistically 
better than other models or the analytic VDFs, espe- 
cially around the peak and the tail. We performed the 
likelihood-ratio test and found that our model fits sig- 
nificantly better for all Rhapsody halos than the SHM 
or the double power-law model at all four radii shown in 

Fig.[U 

In Fig. [2] we also compare three analytic VDFs. For 
the isotropic model shown, the analytic VDF is given 
by Eddington's formula, which gives a one-to-one corre- 
spondence between the density profile and the VDF. For 
anisotropic systems, one must also model the anisotropy 
parameter, defined as fi = 1 — (er 2 + a^)/(2a 2 ), where 

er 2 is the variance in each velocity component. There 
is currently no analytic VDF whose anisotropy profile 
matches that measured in simulations, so we choose three 
simple and representative anisotropic models: constant 
anisotropy (wi th = and 1/2) and t he Osipkov- 
Merritt model (|Osipkovlll979l : lMerritull985l ). The phase- 
space distrib utions of these models ca n be determined 
numerically ([Binnev fc Tremainei |2008|) . For all three 
cases, we adopt the NFW profile as in Eq. ([2]), with the 
best-fit scale radius. For the Osipkov-Merritt model, we 
use the best-fit anisotropy radius. It is shown in Fig. [2] 
and also suggested by the chi-square test for the models 
considered that the analytic VDFs do not describe the 
simulated VDF well. 

Our VDF model, Eq. ([1]), consists of two terms: the 
exponential term and the cut-off term. The origin of the 



the exponential term can be explained by the anisotropy 
in velocity space. Fig. [3] shows the distributions, the dis- 
persion, and the kurtosis of the velocity vectors along 
the three axes of the spherical coordinate. Kurtosis is a 
measure of the peakedness of a distribution, defined as 
(Si v t)/(J2i v i) 2 ~ 3, where Vi is the velocity of the z-th 
particle along one axis, and this value is zero for the nor- 
mal distribution. The ratios of dispersion between the 
three axes are close to one at small radii, and the ratios 
increase with radius. The kurtosis, on the other hand, 
is in general non-zero and decreases with radius. An 
important consequence of the non-zero kurtosis is that 
even if the dispersion along the three axes are similar 
(anisotropy parameter /3 ~ 0), the velocity vectors do 
not follow an isotropic multivariate normal distribution 
in any coordinate system (even after a local coordinate 
transformations). In other words, as long as there exists 
either anisotropy or non-zero kurtosis in a certain coordi- 
nate, the norms of the velocity vectors will not follow the 
Maxwell-Boltzmann distribution. Indeed, Fig. [3] shows 
that in the simulations, one always has non-zero kurto- 
sis and/or anisotropy. Other simulations also indicate 
that the ve locity vectors of dark matter particles have 
anisotropy (lAbel et al.ll201lHSparre fc Hansenll2012[ ) and 
non-zero kurtosis (|Vogelsberger et aLll2009D . We further 
found that if the ratios of dispersion between the three 
axes of a multivariate normal distribution are around 0.2 
to 0.6, the norms of those random vectors will follow a 
distribution which resembles our model without the cut- 
off term, v 2 e xp(— v/vp). (For a for mal discussion on this 
topic, see e.g. lBiornson et ai1l2009l ) This suggests that if 
one can find a coordinate system where the distributions 
of the velocity components are all distributed normally 
(with zero kurtosis), there will be a larger difference be- 
tween the dispersion along the three axes in this new 
coordinate system than in the spherical coordinate. 

The (v 2 sc — v 2 ) p term in our VDF model introduces a 
cut-off at the escape velocity. It further suppresses the 
VDF tail more than the exponential term alone does. De- 
spite that this cut-off term has the form of a power-law 
in (binding) energy, the best-fit values of the parameter 
p does not necessarily reflect the "asymptotic" power- 
law index fc, defined as fc = lim£_>o(^m//<im£), where 
/(£) is the (binding) energy distribution function. The 
relation between fc and the outer density s l ope has been 
studi ed in the literature (jEvans fc Anll2006l : iLisanti et al.l 
2011). However, because dln//dln£ deviates from its 
asymptotic value fc rapidly as £ deviates from zero, 
the asymptotic power-law index fc could be very differ- 
ent from the best-fit power-law index for the VDF tail 
(e.g. v > 0.9w osc ). Furthermore, the shape of the VDF 
power-law tail could be set by recently-ac creted subha- 
jos th at have not been fully phase-mixed (jKuhlen et al.1 
l2012f ) , and hence has no simple relation with the density 
profile. In high-resolution simulated dark matter halos, 
particles stripped off of a still-surviving subhalo are seen 
to significantly impact the tail of the VDF. A larger sam- 
ple of simulations at higher resolution than we consider 
in the current analysis will be needed to further test this 
hypothesis. 

4. HALO-TO-HALO SCATTER IN VELOCITY 
DISTRIBUTIONS 
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Figure 3. Left: The histograms of v r and vj, of the same halo 
shown in Fig. [2] with the best-fit normal distributions (red lines). 
Right: The velocity dispersion cr v /v eB c and the kurtosis, along the 
three axes: v r (red), vg (green), and (blue). Both the disper- 
sion and the kurtosis are measured in spherical shells at different 
r/r s and averaged over all halos in Rhapsody, with the error bars 
showing the 68% halo-to-halo scatter. The dashed lines are only 
to guide the eyes. 



We demonstrated above that there exists a similarity 
in VDFs for a wide range of simulated dark matter halos; 
Eq. (Q~]) provides a good description of this similarity. We 
now quantify explicitly how the VDF depends on r/r s 
and the associated halo-to-halo scatter. Fig. [4] shows 
a scatter plot of the velocity distributions for different 
halos, characterized by the two parameters of Eq. (JT|), for 
different r/r s . The regions of (v ,p) parameter space for 
different r/r s are distinct, which implies that r/r s is the 
most relevant quantity in determining the shape of the 
velocity distribution. We also found that the parameter 
Vo/v esc has a linear relationship in log(r/r s ), as shown in 
the inset of Fig. fJJ 

We note that there is significant degeneracy between 
the two parameters (vq,p). This degeneracy comes from 
the fact that a larger value of p is needed to steepen the 
tail of the VDFs which have larger values of vq. In our 
fitting process we left both parameters free because there 
is no simple relation between vq and p for all radii. Be- 
cause of this degeneracy, there also exists a linear relation 
between p and log(r/V s ). However, since the best-fit p 
is not well-constrained due to the low number of parti- 
cles in the tail of the VDF, the relation between p and 
log(r/V s ) is not well determined cither. 

In Fig. @] we see there exists halo-to-halo scatter even 
for a fixed r/r s . This intrinsic scatter could arise from 
the statistics of the samples or some other physical quan- 
tities. Fig.[S]shows the best-fit vo /v esc at different radii as 
a function of concentration, halo shape (c/a), formation 
time (zi/2), and local density slope (— dlnp/dlnr) re- 
spectively, as defined in lWu et"aLl (|2012f) . We found that 
at a given r/r s (a fixed color), vq/v csc does not have a 
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Figure 4. Distribution of the best-fit parameters, vo and p of 
Eq. JTJ, which describes the simulated VDFs. Each dot represents 
one halo from the RHAPSODY simulation at a certain r/r a : (from 
left to right) 0.15 (blue), 0.3 (red), 0.6 (green), 1.2 (magenta). The 
cross symbols show the best-fit parameters to isotropic analytic 
VDFs obtained from Eddington's formula at corresponding radii. 
The typical uncertainty of the fit is shown in the lower left corner. 
The lower right inset shows the linear relation between vq / v esc and 
log(r/r s ). 
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Figure 5. Scatter plots of the best-fit parameter vo/vesc with 
concentration, halo shape (c/a), formation time (21/2), and local 
density slope on the x-axes respectively. Each dot represents one 
halo from the RHAPSODY simulation at a certain r/r s : (from bot- 
tom to top in each panel) 0.15 (blue), 0.3 (red), 0.6 (green), 1.2 
(magenta). For any fixed r/r s , there is no significant correlation 
between vq/v cbc and the aforementioned quantities on the x-axes. 
See text for details. 



significant correlation with the physical quantities on the 
a;-axis (except for Zi/ 2 in the two smallest radial bins). 
This reinforces the main result of this study: the VDF is 
mostly determined by r/r s (i.e. the gravitational poten- 
tial). We note that the lower left panel of Fig. [S] shows 
a weak correlation between z 1 / 2 and vq/v csc ; however if 
the halos with zi/ 2 < 0.25 are removed, this correlation 
is no longer statistically significant. Halos with recent 
accretion tend to have larger deviations from the NFW 
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profile, and this results in a slight overestimate of the 
best-fit scale radius (fit to an NFW profile). We do not 
expect the Milky Way has had a recent major merger 
with Z1/2 < 0.25. This indicates that for possible Milky 
Way host halos, one can exclude these systems with re- 
cent major mergers, and there will be no remaining cor- 
relation between formation time and vq . 

For Milky Way size halos, it has been suggested 
that the VDF has a universal shape depending only 
on t he velocity disper sion and the local density 
slope ([Hansen et al.lf2006l ). This is related to our finding 
in a way that the magnitude of the velocity dispersion is 
roughly proportional to v csc and the local density slope 
for an NFW profile is given by a monotonic function of 

f J T g 

d\np _ l + 3(r/r.) . 
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However, our study suggests that r/r s is a more funda- 
mental quantity than the local density slope in determin- 
ing the shape of the VDF. Fig. [5] illustrates that vo/v esc 
does not grow with the local density slope when one only 
looks at a fixed r/r s (points with the same color), but it 
does grow with r/r s when the local density slope is fixed. 

5. IMPLICATIONS FOR DIRECT DETECTION 
RATES 

Given the known dependence on r/r s , we can now ex- 
amine the impact on direct dark matter detection ex- 
periments. The differential event rate per unit detector 
mass of dark matter interactions in direct detection ex- 
periments is 



dR 
dQ 



Q 



Poop 
2/j 2 m dr 



-A 2 \F(Q)\< 
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/(V + V e ) 
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where Q is the recoil energy, po is the local dark matter 
density, <7o is the WIMP-nucleus cross section at zero 
momentum transfer, md m is the WIMP mass, p is the 
WIMP-nucleus reduced mass, A is the atomic number of 
the nucleus, |_F(<5)| 2 is the nuclear form factor, w m ; n = 
(QmAr/2/i 2 ) 1//2 for an clastic collision, / is the VDF in 
the Galactic rest frame, a nd v e is the velocity o f Earth 
in the Galactic rest frame (jLewin fc Smith! [19961) . 

With Eq. (j4j one can calculate the event rate given 
VDF and v m [ n . We calculated this rate for each halo 
using the best-fit exponential model of the VDF, for dif- 
ferent w m i n and different r/r s . The results are shown in 
Fig. [5J where we divided the rate by the rate calculated 
from the SHM with conventional parameters vq — 220 
km/s and v esc = 544 km/s for comparison. 

The rate as a function of v m in behaves very differently 
for different r/r s as shown in Fig. [6[ For low values of 
r/r s , the change in detection rates between experiments 
can be much larger than the predictions of the SHM, e.g., 
the r atio between t he rates of CoGeNT ([Aalseth et al.l 
120111 ) to DAMA-I ([Bernabei et al. 2008) is three times 
larger in our model than in the SHM. This clearly mo- 
tivates efforts to better constrain the scale radius of the 
Milky Way: comparing the scatter coming from mea- 
surements of VDF with the intrinsic physical differences 
among halos, the uncertainty on r/r s appears to be the 
dominant contribution to the uncertainty in event rates, 
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Figure 6. Ratio of detection rate predicted by Eq. JT) with pa- 
rameters obtained from Rhapsody, for different r/r s : (from bot- 
tom to top) 0.15 (blue), 0.3 (red), 0.6 (green), to that of the SHM 
with conventional parameters. Vertical dotted lines show v m - ln 
for different detectors: (from left to right) CoGeNT, DAMA-Na, 
XEN ON, CDMS, DAMA-I expressed in (nuc l eus, threshold en- 
ergy) [fA"alseth et al.ll201U;IBernabei et al1 l2008; Apri le et al.H2011l ; 
ICDMS II Collaboration et al.H20iq) , assuming a WIMP mass of 10 
GeV. The error bars show the 68% halo-to-halo scatter, and those 
with wider caps include the scatter in different directions. The x- 
axis is slightly offset for clarity. The lines which connect the data 
points are only to guide the eyes. 

especially for smaller 7J m i n . 

6. DISCUSSION AND CONCLUSION 

When deducing the direct detection event rate from 
cosmological simulations, the primary sources of uncer- 
tainty arise from: (i) finite particle sampling of the VDF, 
(ii) intrinsic scatter from physical processes that affect 
the VDF during the halo formation process (i.e. the halo- 
to-halo scatter), (iii) the quality of the fit and the validity 
of a smooth model, (iv) the observational constraint on 
r/r s for the Milky Way, (v) the variation of the VDFs in 
various directions at a fixed radius, and (vi) the impact 
of baryons. 

An important outcome of our analysis is that at present 
the scatter from (iv) is significantly larger than the cor- 
responding scatter due to each of (i), (ii), and (iii), com- 
bined, by more than two orders of magnitude. This is 
particularly important given that the observational con- 
straint on the sc ale radius suggests the concentration c = 
r vil /r s is 10-20 ([Klvpin et aij|2002t iDeason et al.ll2012fl 
whic h corresponds to r ^/r K ~ 0.15 — 0. 6 dXue et al. | 
[2001 iGnedin et afll2rMlBrown et al.l[20Tfl iBusha et al.l 
1201 If) . Thus, although the distance from the Earth 
to the Galactic cent er is well known ([Ghez et al.l 120081 : 
IGillessen etlll 120091) . we find that the largest current 
theoretical uncertainty on the VDF is the uncertainty in 
r/r s . 

Our determination of the VDF represents an average 
over a spherical shell. In reality, spherical asymmetry 
and substructures will affect the VDF and result in ad- 
ditional scatter along different directions. In the Rhap- 
sody simulations, if we divide the spherical shell into 
several regions while maintaining enough particles (of 
the order 1000) in each analysis region, we find that 
this directional scatter is comparable to the halo-to-halo 
scatter, and that the combined scatter will be 10 — 40% 
larger than only the halo-to-halo scatter, as illustrated 
in Fig. [6] Similar scatter is also seen in the Aquar- 
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IUS Milky Way simulations (|Vogelsberger et al.l l2009f) . 
This directional scatter will grow at larger radii because 
it is a c onsequence of substructures, tidal effects, and 
streams dHelmi et al.1 120031 I Vogelsberger fc Whltel I20ITI: 
IMacieiewski et al.N2011b IPurcell et al.ll2012fl . At present , 
we have no robust way to relate this scatter to direct ob- 
servables, and in practice this directional scatter may be 
the most important uncertainty in determining the direct 
detection rates once other sources have been minimized. 

We have not yet investigated the impact of baryons; we 
expect that adiabatic contraction of dark matter halos 
would raise the velocity but preserve the shape of the 
VDF, so that our model will serve as a useful tool for 
these studies in the future. Baryonic effects in isolated 
halos hav e been studied in the context of da rk matter 
detection (jBruch et al.l [2009: Li ng et alJ [2~01Q K ); however, 
simulating a statistical sample of halos similar to what we 
consider here, with both sufficient resolution and realistic 
baryonic physics is not yet tractable. 

The results presented here highlight the need to sig- 
nificantly improve the determination of the Milky Way 
scale radius. Although the concentr ation is now only 
weakly constrained with present data (jBusha et a.1.11201 it 
iDeason et al.l 12012). improvements will be forthcoming 
with spectroscopy an d astrometry from large scale sur- 
veys ([An et al.l I2012D . Analysis along these lines will 
usher in a new era of complementarity between astro- 
nomical surveys and particle dark matter constraints de- 
duced from terrestrial experiments. 



This work was supported by the U.S. Department of 
Energy under contract number DE-AC02-76SF00515 and 
by a KIPAC Enterprise Grant. Mao is supported by 
a Weiland Family Stanford Graduate Fellowship. We 
thank the Aspen Center for Physics and NSF Grant 
1066293 for hospitality during the workshop "A Theo- 
retical and Experimental Vision for Direct and Indirect 
Dark Matter Detection" . We thank Anatoly Klypin and 
Joel Primack for providing access to the Bolshoi sim- 
ulations, and Peter Behroozi for the halo catalogs for 
both simulations. We thank Tom Abel, Peter Behroozi, 
Michael Busha, Eduardo Rozo, and Li-Cheng Tsai for 
useful discussion. The Rhapsody simulations were run 
using computational resources at SLAC; we gratefully ac- 
knowledge the support of Stuart Marshall, Amedeo Per- 
azzo, and the SLAC computational team. 

REFERENCES 



Aalseth, C. E., Barbeau, P. S., Bowden, N. S., et al. 2011, 

Physical Review Letters, 106, 131301 
Abel, T., Hahn, O., & Kaehler, R. 2011, arXiv: 111 1.3944 
An, J., Evans, N. W., & Deason, A. J. 2012, MNRAS, 420, 2562 
Angloher, G., Bauer, M., Bavykina, I., et al. 2012, European 

Physical Journal C, 72, 1971 
Aprile, E., Arisaka, K., Arneodo, F., et al. 2011, Physical Review 

Letters, 107, 131302 
Behro ozi. P. S.. W echsler. R. H., & Wu, H.-Y. 2011, 

larXiv:1110.43T2l 



Bernabei, R., Belli, P., Cappella, F., et al. 2008, European 

Physical Journal C, 56, 333 
Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second 

Edition, by James Binney and Scott Tremaine. ISBN 

978-0-691-13026-2 (HB). Published by Princeton University 

Press, Princeton, NJ USA, 2008., 
Bhattachariee, P., Chaudhury, S., Kundu, S., & Majumdar, S. 

2012, arXiv:1210.2328 
Bjornson, E., Hammarwall, D., & Ottersten, B. 2009, IEEE 

Transactions on Signal Processing, 57, 4027 

Brown, W. R., Geller, M. J., Kenyon, S. J., & Diaferio, A. 2010, 
AJ, 139, 59 

Bruch, T., Read, J., Baudis, L., & Lake, G. 2009, ApJ, 696, 920 
Busha, M. T., Marshall, P. J., Wechsler, R. H., Klypin, A., & 

Primack. J. 2011, ApJ, 743, 40 
CDMS II Collaboration, Ahmed, Z., Akerib, D. S., et al. 2010, 

Science, 327, 1619 
Deason, A. J., Belokurov, V., Evans, N. W., & An, J. 2012, 

MNRAS, 424, L44 
Eddington, A. S. 1916, MNRAS, 76, 572 
Evans, N. W., & An, J. H. 2006, Phys. Rev. D, 73, 023524 
Feng, J. L. 2010, ARA&A, 48, 495 

Frandsen, M. T., Kahlhoefer, F., McCabe, C, Sarkar, S., & 

Schmidt-Hoberg, K. 2012, J. Cosmology Astropart. Phys., 1, 24 

Ghez, A. M., Salim, S., Weinberg, N. N., et al. 2008, ApJ, 689, 
1044 

Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 
1075 

Gnedin, O. Y., Brown, W. R., Geller, M. J., & Kenyon, S. J. 

2010, ApJ. 720, L108 

Green, A. M. 2012, Modern Physics Letters A, 27, 30004 
Hansen, S. H., Moore, B., Zemp, M., & Stadel, J. 2006, J. 

Cosmology Astropart. Phys., 1, 14 
Helmi, A.. White, S. D. M., & Springel, V. 2003, MNRAS, 339, 

834 

Klypin, A., Zhao, H., & Somerville, R. S. 2002, ApJ, 573, 597 
Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 
102 

Kuhlen, M., Lisanti, M., & Spergel, D. N. 2012, Phys. Rev. D, 86, 
063505 

Kuhlen, M., Weiner, N., Diemand, J., et al. 2010, J. Cosmology 

Astropart. Phys., 2, 30 
Lewin, J. D., & Smith, P. F. 1996, Astroparticle Physics, 6, 87 
Ling, F.-S., Nezri, E., Athanassoula, E., & Teyssier, R. 2010, J. 

Cosmology Astropart. Phys., 2, 12 
Lisanti, M., Strigari, L. E., Wacker, J. G., & Wechsler. R. H. 

2011, Phys. Rev. D, 83, 023519 
Lynden-Bell, D. 1967, MNRAS, 136, 101 

Maciejewski, M., Vogelsberger, M., White, S. D. M., & Springel, 

V. 2011, MNRAS, 415, 2475 
McCabe, C. 2010, Phys. Rev. D, 82, 023530 
Merritt. D. 1985, AJ, 90, 1027 

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 
563 

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 
493 

Navarro, J. F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 
402, 21 

Osipkov, L. P. 1979, Pis ma Astronomicheskii Zhurnal, 5, 77 
Purcell, C. W., Zentner, A. R., & Wang, M.-Y. 2012, J. 

Cosmology Astropart. Phys., 8, 27 
Reed, D. S., Koushiappas, S. M., & Gao, L. 2011, MNRAS, 415, 

3177 

Sparre, M., & Hansen, S. H. 2012, J. Cosmology Astropart. 
Phys., 10, 49 

Stiff, D., & Widrow, L. M. 2003, Physical Review Letters, 90, 
211301 

Strigari, L. E., & Trotta, R. 2009, J. Cosmology Astropart. Phys., 
11, 19 

Vergados, J. D., Hansen, S. H., & Host, O. 2008, Phys. Rev. D, 
77, 023509 

Vogelsberger, M., Helmi, A., Springel, V., et al. 2009, MNRAS, 
395 797 

Vogelsberger, M., & White, S. D. M. 2011, MNRAS, 413, 1419 
Wu, H.-Y., Hahn, O., Wechsler, R. H., Mao, Y.-Y., & Behroozi, 

P. S. 2012, arXiv:1209.3309 
Xue, X. X., Rix, H. W., Zhao, G., et al. 2008, ApJ, 684, 1143 



